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A  pore-network  model  is  developed  to  simulate  liquid  water  transport  in  a  hydrophobic  gas-diffusion 
layer  (GDL)  during  the  operation  of  polymer  electrolyte  membrane  fuel  cells  (PEMFCs).  The  steady  sat¬ 
uration  distribution  in  GDLs  is  determined  through  a  numerical  procedure  using  a  pore-network  model 
combined  with  invasion-percolation  path-finding  and  subsequent  viscous  two-phase  flow  calculation. 
The  simulation  results  indicate  that  liquid  water  transport  in  hydrophobic  GDLs  is  a  strongly  capillary- 
driven  process  that  almost  reaches  the  pure  invasion-percolation  limit  with  zero  capillary  number.  A 
uniform  flux  condition  is  found  to  better  reflect  the  actual  phenomenon  occurring  at  the  inlet  boundary 
for  liquid  water  entering  a  GDL  than  a  uniform  pressure  condition.  The  simulation  further  clarifies  the 
effect  of  the  invaded  pore  fraction  at  a  uniform-flux  inlet  boundary  in  modifying  water  transport  in  GDLs. 
Finally,  the  effect  of  the  GDL  thickness  on  the  steady  saturation  distribution  is  investigated. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  polymer  electrolyte  membrane  fuel  cell  (PEMFC)  is  an  effi¬ 
cient  and  clean  energy-conversion  device  that  generates  electrical 
power  through  the  direct  electrochemical  combination  of  fuel  and 
an  oxidant  [1,2].  Polymer  electrolyte  membrane  fuel  cells  are  oper¬ 
ated  at  relatively  low  temperatures  of  around  80  °C,  which  enables 
good  transient  characteristics,  that  include  short  start-up  times  and 
a  fast  response  to  load  change.  Thus,  PEMFCs  are  considered  to 
be  the  most  suitable  fuel  cell  technologies  for  automotive  power 
generation,  an  application  that  requires  a  high  power  density  and 
good  dynamic  behaviour.  Lately,  water  and  heat  management  has 
attracted  much  research  as  a  way  to  enhance  the  system-level  per¬ 
formance  of  PEMFCs  [3,4]. 

Flooding  in  hydrophobic  gas-diffusion  layers  (GDLs)  is  one  of  the 
most  important  water-management  issues  that  is  closely  connected 
to  the  performance  of  PEMFCs  when  operating  at  high  current  den¬ 
sities.  Accordingly,  many  numerical  studies  have  been  conducted  to 
understand  water  transport  behaviour  and  to  predict  the  saturation 
distribution  in  the  GDLs  of  PEMFCs  [5-14].  It  should  be  noted  that 
these  numerical  models  were  based  on  continuum  transport  the¬ 
ories  for  two-phase  flow  in  porous  media.  Extensive  experimental 
research  has  also  been  performed  to  visualize  liquid  water  trans¬ 
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port  or  to  quantitize  the  liquid  water  content  in  GDLs  by  means 
of  microscopic  observation  [15,16],  magnetic  resonance  imaging 
(MRI),  X-ray  tomography  [17,18],  and  neutron  radiography  [19,20]. 

High-resolution  neutron  radiography  studies  have  been 
reported  very  recently  [21-23],  wherein  the  liquid  water  content 
in  GDLs  was  resolved  at  scales  of  several  micrometers.  The  sat¬ 
uration  distribution  obtained  from  the  resulting  images  [21-23] 
suggested  that  more  complex  processes,  in  addition  to  continuum 
two-phase  flow,  might  operate  during  the  water  transport  in  GDLs. 
Fundamental  modelling  studies  have  also  recently  been  conducted 
to  gain  a  better  knowledge  of  water  transport  in  GDLs  based  on  the 
full-morphological  approach  [24],  the  Lattice-Boltzmann  method 
[25],  and  the  pore-network  model  [26-30]. 

Recent  pore-network  studies  have  shown  that  water  transport  in 
a  hydrophobic  GDL  is  a  strongly  capillary-driven  process  [28-30].  In 
fact,  this  result  is  consistent  with  the  phase  diagram  for  two-phase 
drainage  flow  proposed  by  Lenormand  et  al.  [31  ].  According  to  the 
work  reported  in  [31,32  ],  invasion-percolation  with  capillary  finger¬ 
ing  is  expected  to  be  a  main  transport  mechanism  for  two-phase 
liquid  water  transport  in  GDLs  during  the  operation  of  PEMFCs. 
Invasion-percolation  water  transport  results  in  a  concave  satura¬ 
tion  distribution  along  the  flow  direction  [32],  which  is  different 
from  the  convex  saturation  distribution  predicted  by  conventional 
continuum  two-phase  flow  models. 

In  this  investigation,  a  simplified  pore-network  model  is  devel¬ 
oped  to  simulate  efficiently  steady  water  transport  in  hydrophobic 
GDLs  by  combining  invasion-percolation  path-finding  and  subse- 
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quent  viscous  two-phase  flow  calculation.  The  long  calculation 
time  required  for  full  transient  pore-network  simulations  could  be 
avoided  in  the  simplified  pore-network  model,  which  facilitates 
evaluation  of  the  effects  of  various  parameters  in  shorter  times. 
Utilizing  the  model,  this  study  first  investigates  the  validity  of  two 
inlet  boundary  conditions  for  liquid  water  entering  the  GDLs  as  dis¬ 
cussed  by  Lee  et  al.  [30],  i.e.,  a  uniform  pressure  condition  and  a 
uniform  flux  condition.  Then,  the  effect  of  liquid  water  flux  is  stud¬ 
ied  to  examine  the  dominance  of  capillary  forces  over  viscous  forces 
during  water  transport  in  hydrophobic  GDLs.  The  water  manage¬ 
ment  role  of  microporous  layers  (MPLs),  recently  proposed  by  Nam 
et  al.  [33]  and  Gostick  et  al.  [34],  is  also  validated  by  varying  the 
invaded  pore  fraction  at  a  uniform-flux  inlet  boundary  for  liquid 
water  entering  the  GDLs.  Finally,  the  effect  of  the  GDL  thickness  on 
the  steady  liquid  water  saturation  distribution  is  investigated. 

2.  Theory  and  calculation 

2.1.  Liquid  water  transport  in  GDLs 

Transient  liquid  water  transport  in  a  hydrophobic  GDL  during 
the  operation  of  a  PEMFC  is  a  two-phase  drainage  process  where 
an  invading  non-wetting  fluid  (liquid  water)  displaces  a  defending 
wetting  fluid  (air)  that  initially  filled  the  porous  GDL.  This  process 
is  characterized  by  a  small  capillary  number  (Ca  =  /xinvq/cr  cos  6)  of 
about  10_7-10-8  and  a  moderate  viscosity  ratio  (M  =  /xinv//xdis)  of 
around  10  [28-30].  Such  a  small  capillary  number  is  due  largely 
to  the  small  volume  flux  of  liquid  water,  qinj-,  injected  into  a  GDL 
that,  in  turn,  is  proportional  to  the  water  generation  rate  in  a  cat¬ 
alyst  layer  (CL).  For  example,  a  current  density  of  1 A  cm-2  roughly 
corresponds  to  a  qinj  of  about  10-6  m  s-1 . 

The  small  value  of  Ca  encountered  during  the  operation  of  a 
PEMFC  suggests  that  capillary  forces  rather  than  viscous  forces 
dominate  liquid  water  transport  in  a  GDL.  According  to  the  phase 
diagram  for  the  drainage  processes  [31],  liquid  water  transport  in 
a  hydrophobic  GDL  corresponds  to  the  capillary  fingering  regime, 
wherein  this  regime  is  characterized  by  the  formation  of  long 
and  meandering  transport  paths  due  to  preferential  invasion  of 
non-wetting  fluids  into  the  pores  of  lower  capillary  entry  pres¬ 
sures  [31,32].  Previous  pore-network  studies  [28-30]  have  also 
demonstrated  that  liquid  water  transport  in  a  GDL  is  a  strongly 
capillarity-driven  process  that  is  dominated  by  the  invasion- 
percolation  process  with  capillary  fingering. 

2.1.1.  Capillary  process 

Liquid  water  spontaneously  wets  hydrophilic  capillaries  upon 
contact  due  to  the  adhesive  force  between  the  solid  surface  and 
liquid  water.  The  contact  angle,  0W,  that  liquid  water  makes  with 
air  near  the  solid  surface  is  less  than  90°  for  hydrophilic  capillaries. 
By  contrast,  liquid  water  does  not  spontaneously  wet  hydropho¬ 
bic  capillaries,  and  the  contact  angle  becomes  0W>  90°.  Thus,  the 
liquid  water  pressure,  pw,  is  higher  than  the  ambient  air  pressure, 
pa,  in  hydrophobic  capillaries.  This  also  implies  that  an  additional 
pressure  should  be  applied  to  make  the  liquid  water  penetrate 
hydrophobic  capillaries. 

The  pressure  drop  across  the  water] air  interface  in  circular 
capillaries,  or  the  capillary  pressure,  pc,  is  determined  by  the 
Laplace-Young  equation  [35].  In  this  study,  the  capillary  pressure 
in  hydrophobic  GDLs  is  redefined  for  simplicity  of  calculation  and 
presentation  as: 

Pc  —  Pw  —  Pa  —  cr|COS  6>wl  (  p - ^  ~ET  )  —  - “ -  H) 

V*M  It  2/  G/v/a 

where  a  denotes  the  surface  tension  of  water;  R  i  and  R2  denote  two 
principal  radii  of  curvature;  rw/a  denotes  the  mean  radius  of  curva¬ 


ture  of  the  water|air  interface  defined  as  rw/a  =  2(i?“1  +  R~ 1 )  \  In 
Eq.  (1),  the  capillary  pressure  in  hydrophobic  GDLs  is  treated  as 
positive  and  equal  to  the  gauge  pressure  of  liquid  water  over  the 
ambient  air  pressure.  It  should  also  be  noted  that  the  air  pressure  is 
generally  uniform  in  GDLs  during  the  steady  operation  of  PEMFCs. 

2.1.2.  Condition  of  liquid  water  entering  GDLs 

The  condition  of  liquid  water  at  the  interface  between  a  GDL  and 
a  CL  (or  at  the  interface  between  a  GDL  and  a  MPL)  is  recognized  as 
an  important  factor  for  water  transport  and  saturation  distribution 
in  GDLs  [33,34].  The  following  two  conditions  were  previously  con¬ 
sidered  for  pore-network  studies  of  liquid  water  transport  in  GDLs: 
a  uniform  pressure  condition  [28,29]  and  a  uniform  flux  condition 
[30].  The  former  condition  assumes  that  the  liquid  water  pressure  is 
spatially  uniform  at  the  inlet  boundary  of  a  GDL,  as  if  a  continuous 
liquid  water  reservoir  is  connected  to  all  inlet  pores.  As  a  result,  the 
liquid  water  flux  becomes  spatially  non-uniform  at  the  inlet  bound¬ 
ary  of  the  GDL.  The  uniform  pressure  condition  can  be  implemented 
in  pore-network  models  by  connecting  all  the  inlet  pores  to  a  virtual 
source  pore  that  simulates  the  reservoir. 

On  the  other  hand,  the  uniform  flux  condition  assumes  that  the 
liquid  water  flux  is  spatially  uniform  at  the  inlet  boundary  of  a 
GDL,  as  if  the  liquid  water  is  injected  directly  into  each  inlet  pore. 
Thus,  the  liquid  water  pressure  becomes  spatially  non-uniform  at 
the  inlet  boundary  of  the  GDL.  This  boundary  condition  can  be 
implemented  in  pore-network  models  by  prescribing  a  liquid  water 
injection  rate  for  randomly  selected  inlet  pores.  Lee  et  al.  [30]  dis¬ 
cussed  the  validity  of  the  two  boundary  conditions  and  suggested 
that  a  uniform  flux  condition  is  more  consistent  with  the  uniform 
generation  of  liquid  water  in  a  CL  and  subsequent  invasion  into  a 
GDL.  It  was  also  shown  that  the  uniform  flux  condition  resulted  in 
a  longer  transient  period  toward  a  higher  steady  saturation  level  in 
GDLs  compared  with  the  uniform  pressure  condition. 

2.2.  Pore-network  model 

Liquid  water  transport  in  hydrophobic  GDLs  is  different  from 
that  in  a  bundle  of  straight  hydrophobic  capillaries.  More  com¬ 
plicated  processes  are  expected  to  occur  in  real  GDLs  because  of 
their  complex  pore  structure,  including  irregular  pore  shapes  and 
multiple  pore  connectivities.  Nevertheless  the  fundamental  mech¬ 
anism  of  capillary  transport  in  GDLs  is  believed  to  be  similar  to  that 
explained  with  a  bundle  of  straight  hydrophobic  capillaries. 

The  pore-network  model  is  a  proven  numerical  method  that 
has  been  used  to  investigate  various  transport  phenomena  in 
porous  media  [36-40].  In  the  pore-network  model,  the  complex 
pore  structure  in  porous  media  is  represented  as  a  collection  of 
cross-connected  pores  and  throats.  Although  simplified,  impor¬ 
tant  characteristics  of  the  pore  structure  can  be  preserved  in  the 
pore-network  model  by  a  discrete  choice  of  pore-network  geome¬ 
tries.  Thus,  pore-network  models  have  also  been  used  in  theoretical 
research  for  liquid  water  transport  in  GDLs,  either  to  determine  the 
capillary  pressure  or  relative  permeability  as  a  function  of  satura¬ 
tion  [26,27],  or  to  predict  the  saturation  distribution  [28-30]. 

2.2.1.  Pore-network  generation 

The  pore-networks  for  GDLs  were  constructed  to  be  composed 
of  cross-connected  box-shaped  pores  and  throats  based  on  the 
prescribed  geometrical  parameters  summarized  in  Table  1.  First, 
cubic  cells  with  a  dimension  of  dceli  (=25  p,m)  were  regularly 
arranged  to  fill  the  calculation  domain  (nx  x  ny  x  nz  in  cell  number, 
or  Lx  x  Ly  x  Lz  in  physical  length)  as  shown  in  Fig.  1(a).  Box-shaped 
pores  were  then  generated  to  have  random  dimensions  of  dpx,  dpy, 
and  dpz,  and  positioned  at  the  center  of  each  cubic  cell,  as  illus¬ 
trated  in  Fig.  1(b).  Each  pore  was  connected  to  its  neighboring 
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Table  1 

Parameters  for  present  pore-network  model. 


Parameter 

Explanation 

Value 

Nx  x  Ny  x  Nz 

Cell  number  of  pore-network 

20  x  20  x  10 

Lx  x  Ly  x  Lz 

Size  of  pore-network 

500  p,m  x  500  |jim  x  250  p,m 

f-cell 

Unit  cell  size 

25  |Jim 

ipxi  ipy>  ipz 

Dimension  of  pore  body 

17.5-22.5  |jim 

L a.  itb 

Dimension  of  throat 

5-17.5  pan 

V'p 

Mean  pore  volume9 

9899  pan3  (SD  783  pan3) 

ft 

Mean  throat  radius9 

5.33  p.m  (SD  1.35  p,m) 

e 

Mean  porosity 

0.633  (SD  0.001) 

I< 

Mean  permeability 

3.39  x  10“12  m2  (SD  7.9  x  10“12  m2) 

cr 

Surface  tension  of  water 

0.0645  Nm-1  at70°C 

Pw 

Density  of  water 

978  kg  m-3  at  70  °C 

Viscosity  of  water 

40.38  x  10“5  Pa  sat  70  °C 

Contact  angle  of  water 

120° 

Qinj 

Injection  flux  of  water 

2  A  cm-2  equivalent 

a  Number- mean  values  (not  volume-mean). 


through  the  inlet  boundary  according  to  either  a  uniform  pressure 
or  uniform  flux  condition.  The  injected  volume  flux  of  liquid  water, 
qin j,  at  the  inlet  boundary  is  prescribed  to  be  equivalent  to  a  cur¬ 
rent  density  of  2  A  cm-2.  A  cyclic  pore  connectivity  is  assumed  for 
the  side  planes  of  the  calculation  domain,  as  shown  in  Fig.  1(a), 
which  infinitely  extends  the  calculation  domain  in  the  x-  and  y- 
directions.  The  domain  size  of  the  pore-networks  is  chosen  to  be 
20  x  20  x  10  after  a  grid  dependence  test,  which  corresponds  to  the 
physical  dimension  of  500  p,m  x  500  p>m  x  250  p,m  in  the  x-,  y-,  and 
z-directions. 

It  should  be  noted  that  the  pore-network  geometries  con¬ 
structed  by  the  aforementioned  generation  rules  are  isotropic  in 
nature.  An  anisotropic  pore-network  can  be  obtained  by  differently 
prescribing  the  unit  cell  size,  the  pore  dimension  and  the  throat 
dimension  for  both  in-plane  (x-  andy-direction)  and  through-plane 
(z-direction)  directions. 


six  pores  through  other  box-shaped  throats  having  random  cross- 
sectional  dimensions  of  dta  and  dtb.  The  throats  were  also  aligned 
to  the  center  of  each  cubic  cell  face.  The  ranges  for  the  pore  dimen¬ 
sions  were  0.7dcen  <  dp l=xyz  <  0.9dcen  (17.5  |xm<  dpj=xyz  <  22.5  p,m), 
and  those  for  the  throat  dimensions  were  0.2dcen  <  dti  =  ab  <  0.7dcen 
(5  [xm<dti=ab  <  17.5  p,m).  Similar  pore-network  geometries  have 
also  been  used  in  other  studies  [26,38]. 

The  boundary  conditions  required  for  the  pore-network  calcu¬ 
lation  are  explained  in  Fig.  1(a).  An  inlet  and  outlet  boundary  is 
located  at  the  bottom  and  top  plane  of  a  pore-network,  respec¬ 
tively.  Liquid  water  is  assumed  to  be  introduced  into  a  pore-network 


(a) 

Regularly  stacked 
cubic  cells 


Outlet  (top  plane) 


\ 

Z! _ 

"Periodic 
boundary 
(side  planes) 


Throat  entry 
pressure 

 2c|cos0w| 


Pi  = 


bsi 

1  am  \ 

rp,max  =  ^  min(‘^px  ’  ^py  >  ^pz  ) 


Air  area  Aa 

'  Water  area 
Water 

conductance 


8M. 


Fig.  1.  Pore-network  generation:  (a)  regularly  arranged  cubic  cells,  (b)  box-shaped 
pores  and  throats,  (c)  liquid  water  droplet  in  box-shaped  pore,  and  (d)  invading  front 
of  liquid  water  through  hydrophobic  rectangular  throat. 


2.2.2.  Pore-network  parameters 

Several  geometrical  parameters  required  for  subsequent  calcu¬ 
lation  steps  were  evaluated  just  after  the  pore-network  generation. 
The  pore  volume,  Vp,  of  a  box-shaped  pore  was  calculated  as: 

l/p  =  dpx  X  dpy  X  dp z,  (2) 


and  the  radius  of  the  largest  inscribed  sphere,  rP)max,  was  deter¬ 
mined  as: 

Fp.rnax  =  —  min(dpx,  ^py ’  dpz)-  (3) 

In  Fig.  1(c),  the  inscribed  sphere  radius,  rP)max,  corresponds  to 
the  largest  mean  radius  of  curvature  of  the  water] air  interfaces, 
rw/a,  attainable  in  box-shaped  pores.  Thus,  the  capillary  pressure 
inside  box-shaped  pores  has  the  lowest  value  at  the  saturation  level 
corresponding  to  the  situation  shown  in  Fig.  1(c).  The  transition 
saturation,  Str,  is  thus  defined  as: 


4jZ  rp,min 

Str  =  ^^T' 


(4) 


which  is  the  critical  saturation  for  the  minimum  capillary  pressure 
for  a  pore. 

For  throats,  the  cross-sectional  area,  At,  was  calculated  as: 


At  —  dta  x  dtb,  (5) 

and  the  volume  of  a  throat,  \7t  ij,  connecting  pores  i  and  j  was  deter¬ 
mined  as 


Vt.ij  =  \ij  X  /t.ij ,  (6) 

where  At  ij-  and  /t  ij  are  the  cross-sectional  area  and  the  length  of  the 
throat,  respectively.  From  the  pore-network  generation  rules,  the 
length  of  a  throat,  /t  ij-,  was  calculated  as: 

^t,ij  —  ^cell  —  2^Pi  _  ^pj)»  (7) 

where  dpi  or  dpj  denotes  a  dimension  of  pores  in  the  corresponding 
x-,y-,  or  z-direction. 

The  hydraulic  radius  is  a  parameter  that  influences  both  capillary 
and  viscous  transport  through  rectangular  throats.  The  hydraulic 
radius  of  a  throat,  rt,  is  determined  as: 


2At  _  dtadtb 

Pt  dta  +  dtb 


(8) 


where  Pt  denotes  the  perimeter  of  the  throat,  2(dta  +  dtb).  This 
hydraulic  radius  is  expected  to  be  similar  to  the  mean  radius  of  cur¬ 
vature,  rw/a,  formed  by  liquid  water  invading  a  rectangular  throat, 
as  shown  in  Fig.  1(d). 

The  interface  at  the  invading  non-wetting  phase  front  is  called 
the  main  terminal  meniscus  (MTM)  [36].  As  shown  in  Fig.  1  (d),  there 
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is  another  type  of  water]  air  interface  that  is  formed  at  the  corners 
of  a  long  non-circular  throat.  This  interface  is  called  the  arc  menis¬ 
cus  (AM).  Since  one  principal  radius  of  curvature  becomes  infinite 
as  ->  oo,  the  radius  of  the  arc  meniscus,  ram  (the  same  as  the 
remaining  principal  radius,  R2),  is  equal  to  half  the  mean  radius  of 
curvature  as  ram  =  rw/a/2. 


summarized  to  obtain  the  percolated  clusters  required  for  the  uni¬ 
form  flux  condition. 

Note  that  the  path-finding  procedure  determines  both  the 
essential  clusters  of  invaded  pores  and  throats  due  to  invasion- 
percolation  transport,  as  well  as  the  corresponding  capillary 
pressure  distribution  in  those  percolated  pores. 


2.3.  Pore-network  calculation 


The  present  pore-network  model  calculates  the  capillary  pres¬ 
sure  in  pores  as  a  main  variable,  while  the  capillary  entry  pressure 
and  viscous  resistance  in  throats  are  treated  as  varying  transport 
parameters.  The  first  step  is  to  determine  the  percolated  clus¬ 
ters  of  invaded  pores  and  throats  that  constitute  essential  paths 
for  the  liquid  water  exhaust  through  a  pore-network.  Then,  the 
steady  two-phase  flow  through  the  percolated  clusters  is  solved  in 
order  to  determine  the  capillary  pressure  distribution  required  for 
a  given  liquid  water  flux.  Finally,  the  spatial  saturation  distribution 
in  the  pore-network  is  calculated  by  summing  the  liquid  volume  in 
invaded  pores  and  throats. 


2.3.1.  Invasion-percolation  path-finding 

The  path-finding  procedure  was  implemented  by  referring  to  the 
pure  invasion-percolation  transport  of  liquid  water  in  hydrophobic 
porous  media.  According  to  the  invasion-percolation  rule,  liquid 
water  preferentially  invades  the  throat  with  the  smallest  capillary 
entry  pressure,  pt,  defined  as: 


Pt 


2(7 1  COS  Owl 

n 


(9) 


Note  that  the  hydraulic  radius,  rt,  is  used  to  evaluate  the  mean 
radius  of  curvature  for  the  water] air  interface  of  invading  fronts 
in  rectangular  throats. 

If  certain  pores  and  throats  in  a  pore-network  are  marked 
as  ‘invaded’  at  an  instance  during  the  invasion-percolation  path¬ 
finding  procedure,  then  a  list  can  be  made  of  throats  that  are 
connected  to  the  invaded  pores  but  not  yet  invaded.  This  list  is  called 
the  ‘throat  searching  pool,’  which  registers  all  possible  throat  pas¬ 
sages  available  for  the  next  advance  or  displacement  of  liquid  water. 
Thus,  a  search  for  a  minimum  entry  pressure  should  be  performed 
for  the  throats  in  the  throat  searching  pool. 

Once  a  minimum  capillary  entry  pressure  is  determined,  it  is 
used  to  update  the  capillary  pressure  of  currently  invaded  pores. 
This  pressure  update  is  performed  only  for  the  invaded  pores  whose 
previously  assigned  capillary  pressure  is  smaller  than  the  deter¬ 
mined  minimum  pressure.  Then,  a  throat  corresponding  to  the 
minimum  entry  pressure  and  its  connected  pore  are  marked  as 
invaded  for  the  next  iteration.  The  above  path-finding  procedure 
is  repeated  until  a  throat  that  leads  to  the  outlet  boundary  (the  top 
plane  of  the  domain)  is  finally  invaded.  As  a  final  step,  dead-end 
pores,  which  are  connected  to  only  one  upstream  pore,  are  searched 
in  order  to  check  that  the  capillary  pressures  of  the  dead-end  pores 
are  equal  to  those  of  their  upstream  pores. 

For  the  uniform  pressure  condition,  no  pore  was  marked  as 
invaded  at  the  beginning  of  the  invasion-percolation  path-finding 
procedure.  All  throats  that  were  connected  to  the  inlet  boundary 
(the  bottom  plane  of  the  domain)  were  initially  included  in  the 
throat  searching  pool  and  were  maintained  in  the  pool  until  they 
were  invaded  by  having  a  minimum  entry  pressure  during  the  pro¬ 
cedure. 

For  the  uniform  flux  condition,  a  slightly  different  approach  was 
employed.  A  separate  path-finding  procedure  was  conducted  for 
each  of  the  randomly  selected  inlet  pores  by  marking  a  selected 
pore  as  invaded  (the  throat  between  the  pore  and  the  inlet  bound¬ 
ary  was  initially  also  marked  as  invaded).  Then,  the  invaded  pores 
and  throats  identified  by  the  separate  path-finding  procedures  were 


2.3.2.  Steady  two-phase  flow  calculation 

The  steady  liquid  water  transport  was  calculated  only  for  the 
percolated  clusters  of  invaded  pores  and  throats  determined  by 
the  preceding  path-finding  procedure.  The  viscous  effects  during 
the  two-phase  flow  in  hydrophobic  GDLs  were  fully  considered  by 
including  the  viscous  flow  resistance  of  the  throats.  It  should  be 
noted  that  the  air  pressure,  pa,  was  not  calculated,  but  assumed  to 
be  constant  in  this  study.  This  approach  may  be  justified  from  the 
following  arguments.  Random  stacking  of  fibres  in  GDLs  generally 
results  in  many  small  interstitial  spaces  that  are  difficult  to  be  con¬ 
sidered  in  pore-network  models.  These  spaces  in  hydrophobic  GDLs 
are  too  small  for  liquid  water  to  penetrate,  and  thus,  are  generally 
available  for  air  flow.  In  addition,  the  flow  rate  of  oxygen  during 
the  steady  operation  of  PEMFCs  is  very  small,  leading  to  negligible 
differences  in  the  air  pressure. 

The  conservation  of  volume  at  the  steady-state  is  expressed  as: 

)T  Qm  =  0,  (10) 

j=neiboring  six  pores 


where  Qi^j  is  the  steady  flow  rate  (m3  s-1)  from  a  pore  i  to  its 
neighboring  pore  j.  The  constitutive  equation  for  the  steady  flow 
rate,  Qi_>j,  is  written  as: 

Qi-j  =  Gij[max(pci,  ptiij)  -  max(pcJ,  pUj)],  (11) 


where  pc  i  and  pc  j  are  the  capillary  pressures  in  the  pores,  and  Gy  and 
pt  ij  are  the  flow  conductance  and  capillary  entry  pressure,  respec¬ 
tively,  for  a  throat  connecting  the  i-j  pore  pair.  Eq.  (11 )  reflects  the 
fact  that  liquid  water  transport  occurs  through  a  hydrophobic  throat 
only  when  the  capillary  pressure  is  higher  than  the  entry  pressure 
of  the  throat. 

The  flow  conductance,  Gy,  in  Eq.  (11)  was  derived  from  the 
Hagen-Poiseuille  equation  as  [36]: 


G«  = 


A  -r2 
8/xw/ty 


(12) 


where  Aw  is  the  flow  area  of  liquid  water  in  a  rectangular  throat, 
shown  in  Fig.  1(d),  which  is  calculated  as  Aw  =At-Aa.  The  flow  area 
of  air,  Aa,  was  determined  simply  as: 


=  Him (4  -  Tt). 


(13) 


A  more  complex  expression  for  Aa,  which  considers  the  throat 
geometry  and  the  contact  angle,  can  be  found  in  [36].  Nevertheless, 
Eq.  (13)  was  sufficient  for  the  present  pore-network  study. 

The  radius  of  the  arc  meniscus  in  the  throats,  ram,  can  be  deter¬ 
mined  from  the  capillary  pressure  in  the  upstream  pore,  PcP,  as: 


ram  — 


<7|C0S  Owl 

~W— 


(14) 


A  simpler  approach  was  used  in  this  study,  i.e.,  ram  was  obtained 
by  substituting  p“p  with  the  capillary  entry  pressure  of  the  throat, 
pt.  From  Eq.  (9),  a  constant  ram  is  calculated  for  a  throat  as: 


1 

Jam  =  ^  rt- 


(15) 


Although  Eq.  (15)  neglects  the  variation  of  the  flow  area  in  throats 
according  to  capillary  pressures,  the  changes  in  flow  area  are 
expected  to  be  very  small. 
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In  the  calculation  of  Eq.  (10),  an  alternative  expression  for  Qi^j 
was  used  instead  of  Eq.  (11): 

Qi-j  =  Gy(pcj  - pcj)  +  Gij[max(0,  pt,y  - pCii)  -  max(0, pt>u  - pcJ)]. 

(16) 

The  first  term  in  the  right  hand  side  of  Eq.  (16)  was  included 
in  the  coefficient  matrix  of  the  algebraic  governing  equations, 
while  the  second  term  was  included  in  the  source  vector.  Using 
this  approach,  the  diagonal  dominance  of  the  coefficient  matrix 
required  for  iterative  calculation  was  satisfied  but  a  strong  under¬ 
relaxation  was  needed  due  to  the  highly  non-linear  source  term  in 
Eq-  (16). 

A  constant  volume  flow  rate  of  liquid  water,  Qinj,  was  imposed 
at  the  inlet  boundary  for  calculation  as: 

Qinj  =  ?^(1+2“eff)^d,xy,  (17) 

P\N 

where  Mw  is  the  molecular  mass  of  water  (0.0182  kg  mol-1 );  pw  is 
the  density  of  water,  je  is  the  operation  current  density  (Am-2);  Fis 
the  Faraday  constant  (96,487  C  mol-1 );  aeff  is  the  effective  electro¬ 
osmosis  coefficient;  Adxy  is  the  planar  area  of  the  domain,  Lx  x  Ly. 
For  simplicity  of  presentation,  aeff  was  set  to  0.  For  the  uniform 
flux  boundary,  the  flow  rate  in  Eq.  (17)  was  equally  divided  by  the 
number  of  invaded  inlet  pores,  Ninv,  then  Qi^lN^  was  assigned  to 
each  of  the  invaded  pores.  For  the  uniform  pressure  boundary,  the 
entire  Qinj  was  assumed  to  be  assigned  to  a  virtual  reservoir  cell. 
A  bisect  method  was  used  in  the  calculation  to  find  the  capillary 
pressure  of  the  virtual  reservoir  cell  that  resulted  exactly  in  the 
Outflow  Of  Qinj. 

Once  a  converged  pressure  field  was  obtained,  the  result  was 
used  to  check  if  additional  throats  became  invaded  by  liquid  water. 
If  so,  the  percolated  clusters  were  updated  to  include  those  newly 
invaded  throats  and  also  pores  connected  to  the  throats.  Thus,  the 
capillary  pressure  calculation  is  the  most  computationally  expen¬ 
sive  step  in  the  present  pore-network  model.  On  the  other  hand, 
the  present  pore-network  model  is,  in  fact,  a  simplified  version  in 
comparison  with  previous  transient  pore-network  models  [28-30] 
used  for  water  transport  in  GDLs.  A  fully  transient  model  requires 
a  much  higher  computational  cost  due  to  very  small  time  steps  of 
about  10-6-10-5  s  to  ensure  the  solution  stability,  in  addition  to  a 
relatively  long  physical  time  of  about  60  s  to  reach  the  steady-state 
[30]. 

Transient  phenomena,  such  as  capillary  burst  flow,  cannot  be 
captured  using  the  current  simplified  pore-network  model.  It  is 
considered,  however,  that  the  present  pore-network  model  effi¬ 
ciently  predicts  the  time-averaged  steady  saturation  level  in  GDLs 
with  reduced  computational  costs. 

2.3.3.  Saturation  determination 

At  the  end  of  the  two-phase  flow  calculation,  the  pores  and 
throats  in  a  pore-network  invaded  by  liquid  water  were  identified 
along  with  the  capillary  pressures  within  those  pores.  The  final  step 
was  to  sum  the  liquid  water  volume  in  those  pores  and  throats  to 
obtain  the  saturation  distribution  in  the  pore-network.  For  this,  the 
liquid  water  contents  in  the  cubic  cells  depicted  in  Fig.  1(a)  were 
determined  first.  The  volume  of  liquid  water  in  a  cubic  cell  i,  Wc  it 
can  be  calculated  as: 

Wc,i  =  Wp,j+  Y,  5Wt.ij’  («) 

j=neighboring  six  pores 

where  Wpi  is  the  liquid  water  volume  in  a  pore  i  (the  same  index  is 
used  for  a  cubic  cell  and  a  box-shaped  pore),  and  Wt  y  is  the  liquid 
volume  in  throat-connecting  pores  i  and  j.  A  simple  way  of  deter¬ 
mining  Wp  i  and  Wt  y  is  to  set  Wp ti  =  Vpi  or  Wt  ij  =  Vt  y  when  a  pore  or 
a  throat  is  invaded  by  liquid  water. 


In  this  study,  a  different  approach  was  used  to  determine  Wpi 
and  Wt  y  in  order  to  be  consistent  with  the  pore-network  geome¬ 
tries  shown  in  Fig.  1.  The  liquid  water  volume  in  an  invaded  throat 
was  calculated  as: 


Wt.ij  —  Awjj  x  (19) 

where  Aw  ij  is  the  flow  area  of  liquid  water  in  throats  determined 
during  the  preceding  two-phase  flow  calculation.  The  liquid  water 
volume  in  an  invaded  pore  was  determined  as: 


Wp,i  =  Vp,i  x  Sw(Pc),  (20) 

where  Sw(pc)  is  the  microscale  relationship  between  the  capil¬ 
lary  pressure,  pc,  and  the  saturation,  Sw,  in  the  box-shaped  pores. 
An  approximate  equation  for  pc(Sw)  used  in  Lee  et  al.  [30]  was 
employed  to  obtain  Sw(Pc): 


Pc(3w)  —  2cr|cos  0yj\ 


3_ 

_4tt 


for  Sw  ^  Sf]-. 


(21) 


The  above  equation  was  found  to  slightly  overestimate  Sw  for  a 
given  pc.  Although  a  more  accurate  relationship  could  be  used,  Eq. 
(21 )  was  believed  to  be  sufficient  for  the  present  study.  Note  that  a 
similar  approach  was  also  used  by  Thompson  [39]. 

The  determined  Wc  j  in  the  cubic  cells  was  then  summed  layer- 
by-layer  and  divided  by  the  total  void  volume  in  the  corresponding 
layer.  This  resulted  in  the  layer-by-layer  distribution  of  liquid  water 
saturation.  It  should  be  noted  that  the  current  approach  using 
Eqs.  (19)-(21)  does  not  guarantee  more  accurate  saturation  results 
than  those  obtained  by  the  simpler  method  using  Wp ti  =  Vpi  and 
Wtij  =  Vtij.  In  fact,  the  determined  saturation  generally  has  less 
physical  meaning  in  an  absolute  scale,  but  its  relative  level  is  impor¬ 
tant  when  investigating  the  effects  of  various  conditions  on  water 
transport  and  saturation  distribution. 


3.  Results 

Geometrical  parameters  for  the  randomly  generated  isotropic 
pore-networks  are  summarized  in  Table  1.  The  size  of  the  pore- 
networks  was  20  x  20  x  10,  which  corresponds  to  a  GDL  thickness, 
Lgdl>  of  250  |jim  for  Lcen  =  25  |jim.  The  distributions  of  pore  volume, 
Vp,  throat  area,  At,  transition  saturation,  Str,  and  the  throat  length, 
It,  in  the  simulated  pore-networks  are  presented  in  Fig.  2.  Number- 
mean  values  of  several  parameters  are  also  provided  in  Table  1, 
along  with  their  standard  deviation  (SD).  The  number-mean  throat 
radius,  rt,  is  about  5.33  |xm  in  Table  1 ,  which  corresponds  to  an  area- 
mean  ft  of  5.95  |jim.  These  values  seem  to  be  somewhat  smaller  than 
the  measured  volume-mean  pore  radius  in  GDLs,  rp,  of  about  10  [xm 
[41  ].  Thus,  the  capillary  pressure  is  overestimated  by  about  70%  in 
the  present  pore-network  model  compared  with  actual  pressure  in 
GDLs. 

The  present  pore-network  model  determined  the  mean  sat¬ 
uration  distribution  by  averaging  the  results  of  independent 
simulations  conducted  for  many  different  realizations.  Thus,  a  total 
of  100  realizations  were  simulated  for  each  uniform  flux  or  uniform 
pressure  condition  to  obtain  both  the  mean  saturation  distribution 
and  its  standard  deviation.  Afterwards,  only  10  realizations  were 
simulated  for  each  different  situation  to  determine  only  the  mean 
saturation  distribution  (10  realizations  were  sufficient  to  obtain  a 
converged  saturation  distribution). 

The  mean  porosity,  s,  for  the  simulated  pore-networks  is  esti¬ 
mated  to  be  0.633  (SD  of  0.001 )  and  the  single-phase  permeability, 
I<,  to  be  3.39  x  10-12  (SD  of  0.08  x  10-12)m2.  Williams  et  al.  [41] 
measured  the  through-plane  permeability  of  carbon  paper  (TGP- 
H-120)  with  75.6%  porosity  to  be  about  8.7  x  10-12  m2.  Similarly, 
Feser  et  al.  [42]  reported  that  the  in-plane  permeability  of  carbon 
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Fig.  2.  Statistics  of  pore-network  realizations:  number  frequencies  of  (a)  pore  volume,  Vp;  (b)  throat  area,  At\  (c)  transition  saturation,  Str;  and  (d)  throat  length,  lt. 


paper  (TGP-H-60)  amounted  to  about  5  x  10-12  m2  for  70%  poros¬ 
ity.  Thus,  the  mean  Kof  3.39  x  10-12  m2  seems  to  be  reasonable,  and 
indicates  that  the  present  pore-network  model  properly  considers 
the  pressure  drop  in  GDLs  due  to  viscous  flow  resistance. 

3.1.  Effect  of  inlet  boundary  condition 

The  two  inlet  conditions  describing  the  liquid  water  entering 
hydrophobic  GDLs  were  the  uniform  flux  and  the  uniform  pressure 
conditions.  The  uniform  flux  condition  corresponds  to  a  situation 
where  the  liquid  water  is  introduced  to  a  hydrophobic  GDL  as  if 
it  is  directly  injected  into  each  inlet  pore  by  means  of  many  tiny 
needles.  Thus,  the  liquid  water  flux  is  equally  divided  and  uni¬ 
formly  distributed  to  all  pores  in  the  inlet  boundary.  The  uniform 
pressure  condition  corresponds  to  a  different  situation  where  the 
liquid  water  is  introduced  to  a  hydrophobic  GDL  from  a  continuous 
reservoir  in  contact  with  the  GDL.  By  increasing  the  pressure  of  the 
reservoir,  liquid  water  is  introduced  to  the  hydrophobic  GDL  with 
a  uniform  source-side  pressure  distribution  at  the  inlet  boundary. 

Fig.  3  shows  the  layer-by-layer  distributions  of  liquid  water  sat¬ 
uration  and  capillary  pressure  obtained  by  averaging  the  results  of 
100  realizations  for  each  uniform  flux  or  uniform  pressure  condi¬ 
tion.  The  symbols  in  Fig.  3  denote  the  mean  values  while  the  bared 
ranges  denote  the  confidential  limits  of  ±2  SD  (95.5%).  The  uniform 
flux  condition  results  in  a  higher  saturation  level  and  also  a  higher 
capillary  pressure  than  the  uniform  pressure  condition. 

In  Fig.  3(a),  the  saturation  distribution  with  the  uniform  flux  con¬ 
dition  shows  an  almost  linear  shape  that  connects  the  two  points, 
wherein  Sw  =  0.77  close  to  the  inlet  boundary  and  Sw  =  0.06  close  to 
the  outlet  boundary.  The  average  saturation  for  the  entire  domain 
was  estimated  to  be  0.34.  Note  that  this  saturation  distribution  is 
similar  to  the  steady-state  saturation  distribution  obtained  using 


a  transient  pore-network  model  [30].  On  the  contrary,  the  satura¬ 
tion  distribution  with  the  uniform  pressure  condition  shows  a  more 
concave  shape  in  Fig.  3(a),  where  the  saturation  level  is  estimated  to 
be  0.33  near  the  inlet  boundary  and  0.01  near  the  outlet  boundary. 
Thus,  the  average  saturation  inside  the  entire  domain  is  only  0.10  for 
the  uniform  pressure  condition.  Sinha  and  Wang  [28]  also  predicted 
a  similar  saturation  distribution  based  on  a  transient  pore-network 
model. 

The  capillary  pressure  distribution  also  shows  different  shapes 
for  the  two  boundary  conditions  in  Fig.  3(b).  The  capillary  pressure 
with  a  uniform  pressure  condition  has  an  almost  flat  distribution  at 
around  10  kPa,  except  near  the  outlet  boundary.  On  the  other  hand, 
the  capillary  pressure  distribution  with  a  uniform  flux  condition 
shows  a  gradual  decrease  in  pressure  from  12.3  kPa  near  the  inlet 
boundary  to  10.5  kPa  near  the  outlet  boundary. 

Detailed  transport  behaviour  was  investigated  by  inspect¬ 
ing  the  transport  paths  of  liquid  water  in  a  pore  network.  For 
this  purpose,  an  open-source  Java  viewer,  Jmol,  was  utilized 
(http://www.jmol.org/).  Fig.  4  illustrates  the  percolated  clusters  of 
invaded  pores  and  throats,  which  are  represented  as  balls  and  sticks, 
formed  by  liquid  water  injection  with  the  uniform  flux  condition. 
Similarly,  Fig.  5  shows  the  percolated  clusters  developed  in  a  pore- 
network  according  to  liquid  water  introduction  under  the  uniform 
pressure  condition.  Note  that  two  inlet  conditions  were  applied  to 
an  identical  pore-network  for  direct  comparison.  A  comparison  of 
the  details  given  in  Figs.  4(a)  and  5(a)  indicates  that  a  more  com¬ 
plicated  flow  pattern  develops  in  the  pore-network  with  a  uniform 
flux  condition  than  with  a  uniform  pressure  condition.  Accordingly, 
the  number  of  liquid  water  breakthroughs  is  also  larger  in  Fig.  4(a) 
than  that  in  Fig.  5(a). 

When  liquid  water  is  uniformly  injected  into  each  inlet  pore  of 
a  hydrophobic  GDL  according  to  the  uniform  flux  condition,  it  ran- 
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Fig.  3.  Effect  of  inlet  boundary  condition  on  distribution  of  (a)  liquid  water  satura¬ 
tion  and  (b)  capillary  pressure. 


domly  invades  throats  and  pores  in  search  of  transport  paths  out 
of  the  GDL.  Therefore,  each  of  the  invaded  inlet  pores  has  its  own 
transport  paths  towards  a  liquid  water  breakthrough  point.  Thus, 
a  steady-state  is  achieved  only  when  all  of  the  invaded  inlet  pores 
establish  their  own  steady  transport  paths;  this  results  in  a  longer 
transient  time,  a  higher  saturation  level,  and  a  higher  area-specific 
number  density  of  liquid  water  breakthroughs. 

By  contrast,  the  situation  becomes  very  different  when  liquid 
water  is  introduced  to  a  hydrophobic  GDL  from  a  liquid  water  reser¬ 
voir,  according  to  the  uniform  pressure  condition.  Inlet  pores  with 
larger  dimensions  are  invaded  first  as  the  pressure  of  the  reser¬ 
voir  is  increased,  followed  by  subsequent  random  invasion  of  liquid 
water  through  throats  and  pores.  Once  a  single  breakthrough  path 
is  established,  however,  the  random  invasion  of  liquid  water  stops 
and  all  of  the  liquid  water  flux  is  transported  through  that  break¬ 
through  path.  Thus,  a  steady-state  is  achieved  at  the  instant  that  the 
first  breakthrough  occurs;  this  results  in  a  shorter  transient  time,  a 
lower  saturation  level,  and  a  smaller  area-specific  number  density 
of  liquid  water  breakthroughs. 

In  Fig.  4(a),  the  percolated  clusters  of  invaded  pores  and  throats 
have  been  given  different  colours  according  to  the  corresponding 
breakthrough  points.  Six  breakthrough  paths  identified  in  Fig.  4(a) 
are  individually  plotted  in  Fig.  4(b)  for  better  illustration.  It  can  be 
observed  in  Fig.  4(b)  that  a  breakthrough  point  is  connected  to  sev¬ 
eral  inlet  pores  in  a  certain  region  in  the  inlet  boundary.  If  liquid 
water  is  injected  into  only  one  of  those  inlet  pores,  the  correspond¬ 
ing  breakthrough  point  is  inevitably  invaded.  In  other  words,  the 
injection  of  liquid  water  through  six  deliberately  selected  pores  in 
the  inlet  boundary  can  open  all  six  breakthrough  points  shown  in 


(a) 

Uniform 
flux  BC 


Breakthrough 


Fig.  4.  Liquid  water  transport  paths  with  uniform  flux  condition:  (a)  whole  perco¬ 
lated  clusters  and  (b)  six  individual  clusters. 


Fig.  4(a).  Nevertheless,  the  flow  pattern  in  that  case  is  expected  to 
be  much  simpler  than  that  depicted  in  Fig.  4(a). 

Most  of  the  percolated  clusters  shown  in  Fig.  4(a)  contribute  to 
steady  liquid  water  transport  through  the  pore-network.  By  con¬ 
trast,  most  of  the  percolated  clusters  illustrated  in  Fig.  5(a)  are 
stagnant  (dead-end  clusters),  and  therefore  do  not  contribute  to 
steady  transport.  Fig.  5(b)  presents  only  the  active  clusters  that 
transport  more  than  1%  of  Qinj  introduced  to  the  pore-network, 
which  clearly  indicates  that  only  a  few  percolated  clusters  near  the 
single  breakthrough  point  participate  in  steady  liquid  water  trans¬ 
port.  It  should  also  be  noted  that  the  liquid  water  flux  entering  the 
pore-network  is  highly  non-uniform  in  a  spatial  sense,  wherein  only 
a  small  number  of  throats  seem  to  be  used  for  liquid  water  to  enter 
the  domain. 

The  approximate  regions  of  liquid  water  residing  in  an  identical 
pore-network  are  shown  in  Fig.  6.  A  Boolean  flag  variable,  Finv,  was 
defined  at  the  center  of  each  unit  cell  and  set  to  1  when  the  pore 
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Breakthrough 


(a)  Uniform 
pres.  BC 


(b)  Main  cluesters 

transporting  >1%  Qnj 


Fig.  5.  Liquid  water  transport  paths  with  uniform  pressure  condition:  (a)  whole 
percolated  clusters  and  (b)  active  clusters  transporting  more  than  1%  of  Qmj. 


in  the  unit  cell  was  invaded  (otherwise,  Finv  was  set  to  0).  The  iso¬ 
surfaces  corresponding  to  Finv  =  0.5  were  then  plotted.  Although  the 
presence  of  solid  structures  is  ignored  in  Fig.  6,  these  iso-surfaces 
are  believed  to  illustrate  properly  the  approximate  regions  of  liquid 
water  distribution  in  a  pore-network. 

Fig.  6(a)  corresponds  to  Fig.  4  for  the  uniform  flux  condition,  and 
Fig.  6(b)  to  Fig.  5  for  the  uniform  pressure  condition.  Note  the  iden¬ 
tical  locations  of  liquid  water  breakthrough  points  in  Figs.  6(a)  and 
4,  and  those  in  Figs.  6(b)  and  5.  It  is  demonstrated  again,  in  Fig.  6, 
that  the  uniform  pressure  condition  results  in  a  much  smaller  satu¬ 
ration  distribution  in  the  pore  network  compared  with  the  uniform 
flux  condition.  As  discussed  above,  most  of  the  liquid  water  regions 
shown  in  Fig.  6(a)  contributes  to  the  steady  transport,  whereas  most 
of  the  liquid  water  regions  shown  in  Fig.  6(b)  do  not. 

The  validity  of  the  two  inlet  conditions  has  been  discussed  pre¬ 
viously  by  Lee  et  al.  [30].  It  was  pointed  out  that  a  uniform  liquid 
water  pressure  at  the  inlet  boundary  of  a  GDL  is  difficult  to  attain 
except  when  a  continuous  liquid  water  reservoir  is  in  direct  contact 
with  the  GDL.  Nevertheless,  the  presence  of  a  liquid  water  reser¬ 
voir  just  beneath  a  GDL  is  not  plausible  since  it  will  significantly 
hinder  the  diffusion  of  reactant  gases  towards  a  CL.  In  addition,  the 
spatially  non-uniform  liquid  water  flux  shown  in  Fig.  5(b)  is  incom¬ 
patible  with  the  rather  uniform  generation  of  liquid  water  in  CLs. 
Therefore,  the  uniform  flux  condition  is  believed  to  provide  a  better 
description  of  the  physical  state  of  liquid  water  entering  GDLs  than 
the  uniform  pressure  condition. 

3.2.  Effect  of  liquid  water  flux 

The  distributions  of  liquid  water  saturation  and  capillary  pres¬ 
sure  in  pore-networks  were  obtained  by  varying  the  operation 
current  density,  je,  as  2, 20, 200, 2000,  and  20000  A  cm-2,  as  shown 
in  Fig.  7.  The  operation  current  density  is  directly  proportional  to 
the  liquid  water  flux  through  the  GDL,  as  given  in  Eq.  (17).  For  the 


Fig.  6.  Approximate  regions  of  liquid  water  distribution  with  (a)  a  uniform  flux 
condition  and  (b)  uniform  pressure  condition. 

uniform  flux  condition,  Qmj  in  Eq.  (17)  was  divided  as  Qmj /AW,  and 
then  Qinj/Ninv  was  assigned  to  each  invaded  inlet  pore  that  had  been 
randomly  selected.  Note  that  Fig.  7  was  obtained  with  Ninv  =NxxNy 
(Pinv  =  100%). 

The  saturation  distributions  for  je  of  2  and  20  A  cm-2  are  almost 
identical  in  Fig.  7(a),  and  so  are  the  capillary  pressure  distribu¬ 
tions  in  Fig.  7(b).  These  results  seem  to  be  incorrect  from  the 
viewpoint  of  continuum  two-phase  flow  models.  Recall  that  the 
present  pore-network  model  considers  the  capillary  entry  pres¬ 
sure  and  the  viscous  flow  resistance  as  the  primary  parameters 
governing  liquid  water  transport.  The  liquid  water  flux  influences 
the  viscous  pressure  drop  but  has  no  influence  on  the  capillary 
entry  pressure.  Thus,  the  relative  insensitivity  of  saturation  dis¬ 
tribution  to  liquid  water  flux  can  be  interpreted  as  evidence  that 
the  transport  is  governed  by  capillary  effects  rather  than  viscous 
effects. 

A  current  density  of  2  to  20  A  cm-2  corresponds  to  a  cap¬ 
illary  number,  Ca,  of  about  2.5  x  10_8-2.5  x  10-7.  As  discussed 
earlier,  such  a  small  Ca  with  a  viscosity  ratio,  M,  of  about  10 
corresponds  to  the  capillary  fingering  regime  in  the  phase  dia¬ 
gram  for  two-phase  drainage  flow  [31].  Two-phase  drainage  flow 
in  the  capillary  fingering  regime  is  governed  by  capillary  forces, 
resulting  in  a  concave  saturation  distribution  of  an  invading 
non-wetting  phase.  It  was  proposed  that  the  concave  saturation 
distribution  originates  from  the  invasion-percolation  process  [32]. 
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Fig.  7.  Effect  of  liquid  water  flux,  qinj ,  on  distribution  of  (a)  liquid  water  saturation 
and  (b)  capillary  pressure. 


Fig.  8.  Effect  of  invaded  pore  fraction,  Pin j,  on  distribution  of  (a)  liquid  water  satu¬ 
ration  and  (b)  capillary  pressure. 


At  a  current  density  of  200  A  cm-2,  the  saturation  distribution 
still  maintains  a  rather  concave  shape,  but  its  level  is  slightly 
increased. 

The  saturation  distribution  changes  into  a  convex  shape  when 
the  current  density ,  je,  is  increased  to  2000  and  20000  A  cm-2,  as 
shown  in  Fig.  7(a).  This  transition  in  the  saturation  distribution 
from  a  concave  shape  to  a  convex  shape  indicates  that  the  two- 
phase  drainage  flow  gradually  changes  from  the  capillary  fingering 
regime,  where  capillary  forces  are  dominant,  to  the  stable  displace¬ 
ment  regime,  where  viscous  forces  are  dominant.  The  saturation 
distribution  for  je  of  20000  A  cm-2  has  an  almost  flat  shape  around 
0.8,  which  is  a  characteristic  of  the  stable  displacement  regime.  In 
addition,  the  capillary  pressure  distribution  for  je  of  20000  A  cm-2 
is  also  rather  linear,  which  indicates  the  dominance  of  viscous  flow 
resistance.  Thus,  the  transition  in  the  flow  regime  seems  to  occur 
at  a  je  of  about  200  A  cm-2,  which  corresponds  to  a  Ca  of  about 
2.5  x  10-6.  Sinha  and  Wang  [28]  also  predicted  a  similar  transi¬ 
tion  in  the  saturation  distribution  based  on  the  uniform  pressure 
condition. 

3.3.  Effect  of  invaded  pore  fraction  at  inlet  boundary 

Nam  et  al.  [33]  suggested  that  the  liquid  water  saturation  in 
a  GDL  depends  on  the  number  of  liquid  water  breakthroughs 
at  the  inlet  boundary  of  the  GDL.  Based  on  this  argument,  the 
water  management  role  of  MPLs  was  explained:  MPLs  reduce  the 
area-specific  number  density  of  liquid  water  breakthroughs  from 
CLs  by  merging  many  transport  paths  inside  the  MPLs.  A  simi¬ 
lar  explanation  was  also  recently  proposed  by  Gostick  et  al.  [34] 


after  experimentally  examining  liquid  water  transport  in  GDLs. 
In  order  to  investigate  this  effect,  the  number  of  invaded  inlet 
pores,  Njnj,  was  varied  as  400,  200,  100,  50,  and  25  (for  Ninj-  <400, 
invaded  pores  were  randomly  selected  among  400  inlet  pores). 
The  invaded  pore  fraction  at  the  inlet  boundary,  Pinj,  was  also 
defined  as: 


Thus,  the  above  Ninj-  correspond  to  a  Pinj  of  100,  50,  25,  12.5,  and 
6.25%.  Note  that  Ninj  =400  (Pin j  =  100%)  was  used  for  the  uniform 
flux  condition  in  the  previous  sections. 

Fig.  8(a)  clearly  indicates  that  a  smaller  value  of  Pinj  results  in  a 
reduced  saturation  level  in  GDLs,  especially  near  the  inlet  bound¬ 
ary.  The  saturation  level  near  the  inlet  boundary  was  estimated  to 
be  0.77  for  Pinj  =  100%,  0.68  for  Pinj  =  50%,  0.57  for  Pinj  =  25%,  0.47  for 
Pin j  =  12.5%,  and  0.37  for  Pinj-  =  6.25%.  Pinj  is  observed  to  also  influ¬ 
ence  the  saturation  level  in  far  downstream  regions,  but  the  extent 
of  this  influence  is  not  large.  The  capillary  pressure  distribution  in 
Fig.  9(b)  also  shows  a  similar  trend;  capillary  pressure  near  the  inlet 
boundary  is  more  sensitive  to  Pinj  than  near  the  outlet  boundary.  In 
Fig.  8(a),  the  uniform  pressure  condition  generally  results  in  a  mini¬ 
mum  saturation  level  during  liquid  water  transport  in  hydrophobic 
GDLs. 

Fig.  9  shows  the  percolated  networks  of  invaded  pores  and 
throats  in  a  pore-network  for  a  Pinj  of  100,  25,  and  6.25%.  As  the 
invaded  pore  fraction  increases,  more  complicated  flow  patterns 
appear  to  develop  in  the  pore-network.  This  is  because  more  paths 
are  required  for  liquid  water  transport  according  to  the  increase  in 
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Fig.  9.  Liquid  water  transport  paths  for  invaded  pore  fraction,  Pin j,  of  (a)  100%,  (b) 
25%,  and  (c)  6.25%. 
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the  invaded  pore  number.  The  flow  pattern  near  the  outlet  bound¬ 
ary  is  very  similar,  however,  irrespective  of  the  Piru-  in  Fig.  9.  This 
is  believed  to  be  due  to  the  merging  of  transport  paths  by  the 
invasion-percolation  process  in  the  hydrophobic  pore  network.  In 
fact,  the  number  of  liquid  water  breakthroughs  is  observed  to  be 
less  sensitive  to  P[nj. 

The  results  shown  in  Figs.  8  and  9  support  the  water  manage¬ 
ment  role  of  MPLs  proposed  by  Nam  et  al.  [33  ]  and  Gostick  et  al.  [34]. 
In  addition,  the  dependence  of  liquid  water  saturation  in  GDLs  on 
Pinj  also  explains  the  increased  saturation  level  in  GDLs  during  the 
operation  of  PEMFCs  at  higher  current  densities. 

3.4.  Effect  ofGDL  thickness 

Fig.  10  shows  the  liquid  water  saturation  obtained  by  varying 
the  number  of  unit  cells  in  the  z-direction,  Nz,  as  6,  8,  10,  12,  and 


Fig.  10.  Effect  of  layer  thickness,  Lgd i,  on  distribution  of  liquid  water  saturation  for 
invaded  pore  fraction,  Pinj  ,  of  (a)  100%,  (b)  25%,  and  (c)  6.25%. 


14.  This  pore-network  modification  corresponds  to  the  situation 
where  the  thickness  of  a  GDL,  IG dl,  varies  as  150,  200,  250,  300, 
and  350  p,m.  The  data  in  Fig.  10  clearly  indicate  that  liquid  water 
saturation  depends  on  the  layer  thickness,  wherein  a  smaller  thick¬ 
ness  results  in  a  lower  saturation  level  in  GDLs.  The  saturation  level 
near  the  inlet  boundary  is  less  sensitive  to  the  layer  thickness  com¬ 
pared  with  that  near  the  outlet  boundary.  In  addition,  the  saturation 
distributions  for  Nz  of  12  and  14  are  very  similar.  The  effect  of  the 
invaded  pore  fraction,  Pinj  ,  in  modifying  the  saturation  distribution 
is  also  observed  in  Fig.  10. 

The  continuum  two-phase  flow  theories  generally  correlate  the 
increased  layer  thickness  of  porous  media  to  the  increased  viscous 
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Fig.  11.  Liquid  water  transport  paths  for  layer  thickness,  Lgd i,  of  (a)  150  p,m,  (b) 
250  pan,  and  (c)  350  |jim. 


flow  resistance.  Thus,  a  larger  layer  thickness  requires  a  higher  vis¬ 
cous  pressure  drop,  thereby  resulting  in  a  higher  saturation  level. 
An  alternative  explanation  can  be  made  based  on  the  invasion- 
percolation  mechanism.  The  invasion-percolation  of  liquid  water 
that  searches  for  transport  paths  towards  breakthrough  points  in 
the  outlet  boundary  can  be  minimized  when  the  porous  media 
are  thin.  Accordingly,  the  flow  structure  becomes  simpler  as  the 
thickness  is  reduced,  resulting  in  a  lower  saturation  level  in  GDLs. 

The  percolated  networks  of  invaded  pores  and  throats  for  Nz  of  6, 
10,  and  14  are  shown  in  Fig.  11,  where  the  number  of  breakthrough 
points  increases  as  the  layer  thickness  decreases.  This  result  indi¬ 
cates  that  the  merging  of  transport  paths  appears  to  be  hindered 
in  thin  porous  media.  A  similar  trend  is  also  found  in  Fig.  10  that 
the  saturation  level  near  the  outlet  boundary  increases  as  the  layer 
thickness  decreases.  The  saturation  levels  near  the  outlet  are  esti¬ 


mated  to  be  0.02  for  350  p,m,  0.04  for  300  pan,  0.06  for  250  pm,  0.07 
for  200  pm,  and  0.11  for  150  pm. 

4.  Discussion 

The  present  results  reconfirm  the  findings  of  previous  pore- 
network  studies  [28,30]  that  liquid  water  transport  in  ‘uniformly’ 
hydrophobic  GDLs  is  a  strongly  capillary-driven  process.  As  dis¬ 
cussed  in  [28,30],  this  conclusion  is  also  consistent  with  the  phase 
diagram  for  the  two-phase  drainage  flow  in  porous  media  proposed 
by  Lenormand  et  al.  [31].  The  invasion-percolation  process  [32]  is 
thus  believed  to  be  the  main  transport  mechanism  for  the  capil¬ 
lary  transport  of  liquid  water  in  hydrophobic  GDLs.  In  addition,  the 
uniform  flux  condition  is  also  found  to  provide  a  better  descrip¬ 
tion  of  the  actual  process  occurring  at  the  inlet  boundary  of  GDLs 
compared  with  the  uniform  pressure  condition. 

Accepting  the  dominance  of  capillary  processes,  several  aspects 
of  liquid  water  transport  in  uniformly  hydrophobic  GDLs  can  be 
explained.  A  water  generation  rate  equivalent  to  a  current  density 
of  2-20  A  cm-2  is  found  to  make  almost  no  difference  in  the  steady 
saturation  distribution  in  Fig.  7.  This  result  suggests  that  the  liq¬ 
uid  water  flux  in  the  range  of  normal  PEMFC  operation  has  less 
impact  on  the  saturation  level  in  GDLs.  Instead,  the  data  in  Fig.  8 
point  out  that  the  saturation  level  in  GDLs  depends  on  the  num¬ 
ber  of  invaded  pores  at  the  uniform-flux  inlet  boundary.  The  higher 
saturation  levels  in  GDLs  experimentally  observed  during  the  high 
current  density  operation  of  PEMFCs  can  then  be  explained  by  the 
increased  number  of  liquid  water  breakthroughs  from  CLs. 

The  water  management  role  of  MPLs  recently  proposed  by  Nam 
et  al.  [33  ]  and  Gostick  et  al.  [34]  can  also  be  explained  by  the  depen¬ 
dence  of  the  saturation  level  in  GDLs  on  the  invaded  pore  fraction. 
Many  transport  paths  formed  by  the  relatively  uniform  injection  of 
liquid  water  from  CLs  are  expected  to  merge  inside  the  MPLs  and 
thereby  result  in  a  significant  reduction  in  the  number  of  liquid 
water  breakthroughs  towards  the  GDLs.  This  process  in  MPLs  is  also 
inferred  by  Fig.  4,  where  the  number  of  breakthrough  outlet  pores 
(six)  is  much  smaller  than  that  of  invaded  inlet  pores  (20  x  20). 
Thus,  MPLs  are  believed  to  minimize  the  invaded  pore  fraction  at 
the  inlet  boundary  of  GDLs,  thereby  leading  to  a  reduced  saturation 
level  in  the  GDLs. 

A  smaller  thickness  is  found  to  reduce  the  saturation  level  in 
GDLs,  as  shown  in  Fig.  10,  which  is  understandable  from  the  view¬ 
point  of  invasion-percolation  transport.  When  the  thickness  of  a 
GDL  is  small,  it  reduces  the  invasion-percolation  process  of  liq¬ 
uid  water  that  searches  for  transport  paths  toward  breakthrough 
points.  This  trend  is  also  observed  by  less  complicated  flow  patterns 
obtained  for  smaller  layer  thicknesses,  as  shown  in  Fig.  11 .  The  num¬ 
ber  of  liquid  water  breakthroughs  increases  as  the  thickness  of  GDLs 
decreases  in  Fig.  11 ,  which  is  consistent  with  the  inverse  proportion¬ 
ality  between  the  area-specific  breakthrough  number  density  and 
the  layer  thickness  suggested  by  Nam  et  al.  [33]. 

Although  not  appropriate  as  an  inlet  condition,  the  uniform 
pressure  condition  is  worth  investigating  because  it  leads  to  a  min¬ 
imum  saturation  level  in  GDLs.  The  condensation  of  water  vapour 
inside  GDLs  is  also  important  since  it  generally  results  in  a  higher 
saturation  level  than  that  required  only  for  liquid  water  trans¬ 
port.  Therefore,  further  studies  are  currently  needed  to  resolve 
these  issues  from  the  viewpoint  of  invasion-percolation  transport 
in  GDLs.  In  addition,  the  effect  of  the  mixed  wettability  in  GDLs  due 
to  imperfect  hydrophobic  coating  [29]  also  needs  to  be  clarified  in 
further  studies. 

5.  Conclusions 

Liquid  water  transport  and  saturation  distribution  in  hydropho¬ 
bic  GDLs  has  been  investigated  by  using  a  simplified  pore-network 
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model  combined  with  invasion-percolation  path-finding  and  sub¬ 
sequent  steady  two-phase  flow  calculation.  The  pore-network 
results  reconfirmed  the  dominance  of  the  capillary  process  dur¬ 
ing  liquid  water  transport  in  hydrophobic  GDLs.  The  results  also 
show  that  the  saturation  level  in  GDLs  can  be  lowered  by  reduc¬ 
ing  the  invaded  pore  fraction  or  by  reducing  the  layer  thickness. 
Several  aspects  of  liquid  water  transport  in  uniformly  hydropho¬ 
bic  GDLs  are  explained  in  terms  of  the  invasion-percolation 
mechanism. 
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